home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / demos / ode / ode.dem < prev   
Text File  |  1999-09-16  |  3KB  |  102 lines
  1. mode(5)
  2. //1- Simple Example
  3.  xbasc();
  4.  deff('[xd]=lin(t,x,a)','xd=a*x')
  5.  a=[1 1;0 2];
  6.  ea=ode(eye(2,2),0,1,list(lin,a)),exp(a)
  7.  
  8.  t=0:0.1:3;
  9.  ee=ode(1,0,t,list(lin,1));plot2d1("onn",t',ee',-(1:2),"121",'x1@x2');
  10.  xtitle('dx=a*x','t',' ')
  11.  halt();xbasc();
  12.  
  13. //2- chemical process (stiff)
  14. mode(-1)
  15. titlepage(['Integration of ODE:';...
  16.         'dy1/dt=-0.04*y1 + 1d4*y2*y3';...
  17.         'dy3/dt= 3d7*y2*y2';...
  18.         'dy2/dt= -dy1/dt - dy3/dt';...
  19.         'finding points such that';...
  20.         'y1=1.e-4 or y3=1.e-2'])
  21.  
  22.  rect=[1.d-5,-0.1,1d11,1.1];
  23.  deff('[yd]=chem(t,y)',[
  24.                          'yd(1)=-0.04*y(1) + 1d4*y(2)*y(3);';
  25.                          'yd(3)= 3d7*y(2)*y(2);';
  26.                          'yd(2)= -yd(1) - yd(3);'])
  27.  
  28.  comp(chem)
  29.  mode(1)
  30.  t=[1.d-5:0.02:.4 0.41:.1:4 40 400 4000 40000 4d5 4d6 4d7 4d8 4d9 4d10];
  31.  rtol=1.d-4;atol=[1.d-6;1.d-10;1.d-6];
  32.  y=ode([1;0;0],0,t,rtol,atol,chem);
  33.  halt();xbasc();
  34.  plot2d1("oln",t',(diag([1 10000 1])*y)',-(1:3),"111",' y1@10000 y2@y3',rect)
  35.  halt();
  36.  nt=prod(size(t));
  37.  deff('[y]=surf(t,x)','y=[x(1)-1.e-4;x(3)-1.e-2]')
  38.  comp(surf)
  39.  [y,rd,w,iw]=ode('root',[1;0;0],0,t,rtol,atol,chem,2,surf);rd
  40.  while rd<>[] then ..
  41.   [nw,ny]=size(y);
  42.   k=find(rd(1)>t(1:nt-1)&rd(1)<t(2:nt));..
  43.   write(%io(2),[rd(1);y(:,ny)]','(''t='',e10.3,'' y='',3(e10.3,'',''))')
  44.   plot2d1("oln",rd(1)',(diag([1 10000 1])*y(:,ny))',[3,3,3],"000");
  45.   //mplot([rd(1);rd(1);rd(1)],diag([1 10000 1])*y(:,ny));..
  46.   [y,rd,w,iw]=ode('root',[1;0;0],rd(1),t(k+1:nt),rtol,atol,chem,2,surf,w,iw);..
  47.  end
  48. halt();xbasc();
  49. //implicit
  50. mode(-1)
  51. //
  52. titlepage(['Implicit ODE:';...
  53.           ' ';...
  54.         'dy1/dt=-0.04*y1 + 1d4*y2*y3';...
  55.         'dy2/dt=0.04*y1 - 1d4*y2*y3-3d7*y2*y2';...
  56.         '  1   = y1+y2+y3'])
  57.  
  58.  deff('[r]=chemres(t,y,yd)',[
  59.                          'r(1)=-0.04*y(1)+1d4*y(2)*y(3)-yd(1);';
  60.                          'r(2)=0.04*y(1)-1d4*y(2)*y(3)-3d7*y(2)*y(2)-yd(2);'
  61.                          'r(3)=y(1)+y(2)+y(3)-1;'])
  62.  
  63.  
  64.  deff('[p]=chemad(t,y,p)','p=p+diag([1 1 0])')
  65.  
  66.  deff('[p]=chemjac(t,y,yd)',['p=[-0.04, 1.d4*y(3)        ,1.d4*y(2);';
  67.                             '    0.04,-1d4*y(3)-6d7*y(2),-1d4*y(2);';
  68.                             '    1 ,    1,                  1       ]'])
  69.  
  70. write(6,' macro g(t,y)-a(t,y)*ydot')
  71. write(6,[' ';' macro  x=a(t,y)+p'])
  72. write(6,[' ';'jacobian'])
  73. comp(chemres);comp(chemad);comp(chemjac)
  74. mode(1)
  75. y0=[1;0;0];
  76. yd0=[-0.04;0.04;0];
  77. t=[1.d-5:0.02:.4 0.41:.1:4 40 400 4000 40000 4d5 4d6 4d7 4d8 4d9 4d10];
  78. rtol = 1d-4;atol=[1.d-6;1.d-10;1.d-6];
  79. y=impl(y0,yd0,0,t,rtol,atol,chemres,chemad,chemjac);
  80.  halt();xbasc();
  81. plot2d1("oln",t',(diag([1 10000 1])*y)',-(1:3),"111",'y1@10000 y2@y3',rect)
  82.  
  83. mode(-1)
  84.  halt();xbasc();
  85.  
  86. titlepage('lorenz ode ');
  87.  
  88. deff('[ydot]=lorenz(t,y)',...
  89. "x=y(1);...
  90. a=[-10,10,0;28,-1,-x;0,x,-8/3];...
  91. ydot=a*y")
  92. deff('[j]=jacobian(t,y)',...
  93. "x=y(1);yy=y(2);z=y(3);...
  94. j=[-10,10,0;28-z,-1,-x;-yy,x,-8/3]")
  95. comp(lorenz);comp(jacobian);
  96. y0=[-3;-6;12];t0=0;step=0.01;t1=10;
  97. instants=t0:step:t1;
  98. y=ode(y0,t0,instants,lorenz,jacobian);
  99. xbasc(0);param3d(y(1,:),y(2,:),y(3,:))
  100.  
  101.  
  102.